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Abstract 

Gaussian beams are asymptotically valid high frequency solutions to hyperbolic partial 
differential equations, concentrated on a single curve through the physical domain. They 
can also be extended to some dispersive wave equations, such as the Schrodinger equation. 
Superpositions of Gaussian beams provide a powerful tool to generate more general high 
frequency solutions that are not necessarily concentrated on a single curve. This work is 
concerned with the accuracy of Gaussian beam superpositions in terms of the wavelength 
e. We present a systematic construction of Gaussian beam superpositions for all strictly 
hyperbolic and Schrodinger equations subject to highly oscillatory initial data of the form 
Ae 1 ®^ . Through a careful estimate of an oscillatory integral operator, we prove that the fc-th 
order Gaussian beam superposition converges to the original wave field at a rate proportional 
to e fc//2 in the appropriate norm dictated by the well-posedness estimate. In particular, we 
prove that the Gaussian beam superposition converges at this rate for the acoustic wave 
equation in the standard, e-scaled, energy norm and for the Schrodinger equation in the L 2 
norm. The obtained results are valid for any number of spatial dimensions and are unaffected 
by the presence of caustics. We present a numerical study of convergence for the constant 
coefficient acoustic wave equation in R 2 to analyze the sharpness of the theoretical results. 



1 Introduction 

In simulations of high frequency wave propagation, a large number of grid points is needed to resolve 
and maintain an accurate in time representation of the wave field. Consequently, in this regime, 
direct numerical simulations are computationally expensive and at sufficiently high frequencies, 
such simulations are no longer feasible. To circumvent this difficulty, approximate high frequency 
asymptotically valid methods are often used. One such popular approach is geometrical optics 
[7j [33] , which is obtained in the limit when the frequency tends to infinity. This method is also 
known as the WKB method or ray-tracing. The solution of the partial differential equation (PDE) 
is assumed to be of the form 

a(t,y,e)e^^, (1) 

where 1/e is the large high frequency parameter, 4> is the phase, and a is the amplitude of the 
solution having the Debye expansion in terms of e, a(t,y,e) — ^2j =Q £-'cij(t,y). The phase and 
amplitudes dj are independent of the frequency and vary on a much coarser scale than the full wave 
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solution. They can therefore be computed at a computational cost independent of the frequency. 
However, the geometrical optics approximation breaks down at caustics, where rays concentrate 
and the predicted amplitude is unbounded [2H US]- The consideration of difficulties caused by 
caustics, beginning with Keller in [TB] and Maslov and Fedoriuk (see [25]), led to the development 
of the theory of Fourier integral operators, e.g., as given by Hormander in |10) . 

Gaussian beams form another high frequency asymptotic model which is closely related to 
geometrical optics. However, unlike geometrical optics, Gaussian beams do not breakdown at 
caustics. For Gaussian beams, the solution is also assumed to be of the geometrical optics form 
([1]) , but a Gaussian beam is a localized solution that concentrates near a single ray of geometrical 
optics in space-time. Although the phase function is real-valued along the central ray, Gaussian 
beams have a complex-valued phase function off their central ray. The imaginary part of the phase 
is chosen such that the solution decays exponentially away from the central ray, maintaining a 
Gaussian-shaped profile. To form a Gaussian beam solution, we first pick a ray and solve a system 
of ordinary differential equations (ODEs) along it to find the values of the phase, its first and 
second order derivatives and the amplitude on the ray. To define the phase and amplitude away 
from this ray to all of space-time, we extend them using a Taylor expansion. Heuristically speaking, 
along each ray we propagate information about the phase and amplitude and their derivatives that 
allows us to reconstruct the wave field locally in a Gaussian envelope. The existence of Gaussian 
beam solutions has been known since sometime in the 1960's, first in connection with lasers, see 
Babic and Buldyrev [2J. Later, they were used in the analysis of propagation of singularities in 
partial differential equations by Hormander [TT] and Ralston [30] . 

In this article, we are interested in the accuracy of Gaussian beam solutions to m-th order 
linear, strictly hyperbolic PDEs with highly oscillatory initial data of the type 

Pu = 0, (t, y) G (0, T] x M", (2) 

N 

d e t u(0,y) = e^Xi^'M^ > i = Q,...,m-l, 

3=0 

where the strictly hyperbolic operator, P, is defined in Section [2~T1 the real valued phase, $(y) 
belongs to C°°(Kq; K) for some compact set Kq C 1", and the complex valued amplitudes, Aij(y) 
belong to (Kq;C). Furthermore, we will assume that |V$(y)| is bounded away from zero on 
Kq. As a special case, we include the acoustic wave equation, 

um - c{y fAu = 0, (t, y) G (0, T] x M" , 

N N 

u(0,y) =J2^MM^ {v)/£ , and u t (0,y) = A^{ y y Hv)/e ■ (3) 

3=0 j=0 

We also treat the dispersive Schrodinger equation, 

-i SUt -L.Au + V{y)u = 0, (t, y) G (0, T] x R™, (4) 

N 

u(Q 1 y) = Y j e ] A 3 { y y Hv)/e ■ 

3=0 

As above, we will assume that $ G C°°(Kq;R) and Aj G (Kq\C). Furthermore, we assume 
the potential V(y) is smooth and bounded along with all its derivatives, d@V G Cg°(M. n ) for all j3. 
For the Schrodinger equation, the asymptotic parameter e appears in the equation and there is no 
need to assume the bound on |V$| that is necessary for strictly hyperbolic equations. 
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Since these partial differential equations are linear, it is a natural extension to consider sums 
of Gaussian beams to represent more general high frequency solutions that are not necessarily 
concentrated on a single ray. This idea was first introduced by Babic and Pankratova in [3] and 
was later proposed as a method for wave propagation by Popov in |28j . The sum, or rather the 
integral superposition, of Gaussian beams in the simplest first order form can be written as 



where Kq is a compact subset of R™ and the phase that defines the Gaussian beam is given by 



The real vector p(t; z) is the direction of wave propagation and the matrix M(t; z) has a posi- 
tive definite imaginary part and it gives Gaussian beams their profile. Extensions of the above 
superposition are possible in several directions, including using higher order Gaussian beams in 
the superposition and using a sum of several superpositions to approximate the different modes 
of wave propagation. Higher order Gaussian beams are created by using an asymptotic series for 
the amplitude and using higher order Taylor expansions to define the phase and the amplitude 
functions, see (|16l) and (|17p . Also, for higher order beams, a cutoff function (ITBl is necessary to 
avoid spurious growth away from the central ray. Superpositions with higher order Gaussian beams 
have an improved asymptotic convergence rate. For m-th order strictly hyperbolic PDEs, which 
have m pieces of initial data, we use m different Gaussian beam superpositions chosen in such a 
way so that their sum approximates the initial data. Each of these superpositions corresponds to 
one of the m distinct modes of wave propagation. 

Accuracy studies for a Gaussian beam solution mqb have traditionally focused on how well it 
asymptotically satisfies the PDE, i.e. the size of the norm of Pmgb in terms of e. The question 
of determining the error of the Gaussian beam superposition compared to the exact solution was 
thought to be a rather difficult problem decades ago, see the conclusion section of the review article 
by Babic and Popov [I]. However, some progress on estimates of the error has been made in the 
past few years. This accuracy study was initiated by Tanushev in [34], where a convergence rate 
was obtained for the initial data. Some earlier results on this were also established by Klimes in 
[T8] . The part of the error that is due to the Taylor expansion off the central ray was considered by 
Motamed and Runborg in [27] for the Helmholtz equation. Liu and Ralston [22, 23J gave rigorous 
convergence rates in terms of e for the acoustic wave equation in the scaled energy norm and for 
the Schrodinger equation in the L 2 norm. However, the error estimates they obtained depend on 
the number of space dimensions in the presence of caustics, since the projected Hamiltonian flow to 
physical space becomes singular at caustics. The superpositions in §5§ can also be carried out over 
both x and p in full phase space through the Hamiltonian map, (z,pq) — > (x(t; z,po),p(t; z,po)), as 
shown in [22J 123] ■ I n this formulation the Hamiltonian flow is regular and there are no caustics, so 
we expect to obtain a dimcnsionally independent error estimate. This has been confirmed for the 
wave equation by Bougacha, Akian and Alexandre in [5] for the case of initial data based on the 
Fourier-Bros-Iaglonitzer (FBI) transform, a result which is inspired by the work of Rousse and 
Swart on the Herman-Kluk propagator for the Schrodinger equation |31H32) . From a computational 
stand point, the full phase space formulation is more expensive. 

Building upon these recent advances, together with an application of a non-squeezing argument 
proved in Lemma 02 we are able to provide a definite answer to the question of accuracy for 
Gaussian beam superposition solutions. More precisely, we obtain dimensionally independent 
estimates for the superposition in physical space for general m-th order strictly hyperbolic PDEs 
and the Schrodinger equation. Our main result is the following theorem. 





<t>(t,y,z) = 0o{t;z) +yp(t;z) + y -M(t;z)y . 



(6) 
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Theorem 1 Let u be the exact solution to the PDEs considered, @) ; |3[) and under the stated 
assumptions on initial data and partial differential operators. Moreover, let u k be the corresponding 
k-th order Gaussian beam superposition given in Section and Section \2.Sl with a sufficiently 
small cutoff parameter rj when k > 1. We then have the following estimates. For the m-th order 
strictly hyperbolic PDE |]|[), 

m— 1 
1=0 

For the acoustic wave equation 

\\u(t,-)-u k (t,-)\\ E <C(T)e k ^ , 
where \ \ ■ \\e is the scaled energy norm h21)) . Finally, for the Schrddinger equation |^J), 

\\u(t,-)-u k (t,-)\\ L 2 <C{T)e k / 2 . 

This improves on the results in [22] [23] where the last two error estimates were also proved, but 
with an additional factor e~ 7 in the right hand side, where 7 = (n — l)/4 for the wave equation 
and 7 = n/4 for the Schrddinger equation. Note that the rescaling by e" 1-1 in the first estimate is 
convenient here, since it exactly balances the rate at which the corresponding norm of the initial 
data for the PDE goes to infinity as e — > 0. 

At present there is considerable interest in using numerical methods based on superpositions 
of beams to resolve high frequency waves near caustics, which began in the 1980's with numerical 
methods for wave propagation in [28[ I15[ [6] and more specifically in geophysical applications 
in [17l E E]- Recent work in this direction includes simulations of gravity waves [36], of the 
semiclassical Schrddinger equation [T31 [2D], and of acoustic wave equations [3S[ [M]. Numerical 
techniques based on both Lagrangian and Eulerian formulations of the problem have been devised 
[14l [2T1 [20l [26] . A numerical approach for general high frequency initial data closely related to 
the FBI transform, but avoiding the cost of superposing over all of phase space, is presented in 
|29| for the Schrddinger equation. Numerical approaches for treating general high frequency initial 
data for superposition over physical space were considered in [351 [T] for the wave equation. Our 
theoretical results show that the numerical solutions found in these papers will be accurate when 
£ < 1. 

To test the sharpness of the theoretical convergence rates, we present a short numerical study in 
the case of the acoustic wave equation with constant sound speed. Our numerical results indicate 
that the theoretical rates are sharp for even order k, but similar to the result in [27] . we observe a 
gain in the convergence rate of a factor of e 1 / 2 for beams of odd order k, which suggests that the 
actual convergence rate is ©(e^/ 2 ! ). 

This paper is organized as follows: Section [2] introduces Gaussian beams and their superposi- 
tions for m-th order strictly hyperbolic equations. Furthermore, we construct Gaussian beams for 
the Schrddinger equation. Section [3] is devoted to error estimates for Gaussian beam superposi- 
tions. Detailed norm estimates of the oscillatory operators used in obtaining the error estimates 
are given in Section [4] Numerical validation of our results is finally presented in Section [5] 

2 Construction of Gaussian Beams 

In this section, we outline the construction of Gaussian beam superpositions for strictly hyperbolic 
PDEs. We also construct Gaussian beams for the Schrddinger equation. 
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2.1 Hyperbolic Equations 

Let P = P m + L be a linear strictly hyperbolic m-th order partial differential operator (PDO) in 
n dimensions with 

m-l / \ 

p TO = 9 t m +£ Yl 9p(t,v)ae\d>, (?) 

and L a differential operator of order m — 1. The principal symbol of P, denoted by er m (t, y, r,p), 
is defined by the formal relationship P m = <j m (t,y, —idt, —id y ). Following [13l [30], we make the 
assumptions: 

(HI) The coefficients gp{t,y) are smooth functions, bounded in t and y along with all their deriva- 
tives, dfd£gi3 G C£°(R n ) for all l,a. 

(H2) For \p\ ^ 0, the principal symbol o~ m (t, y, T,p) has m distinct real roots, when it is considered 
as a polynomial in r. 

(H3) These roots are uniformly simple in the sense that 

da m (t,y,r,p) 



> co\p\ m whenever a m (t, y, r,p) = 0. (8) 



We consider a null bicharacteristic (t(s),x(s),r{s),p(s)) associated with the principal symbol 
er m , defined by the Hamiltonian system of ODEs: 

• _ da m . _ da m , _ da m . _ da m 

(h ' X ~~dp~ : T ~~~dT' P ~~~dy~ ' [ ) 

and initial conditions (t(0),x(0),r(0),p(0)) such that cr m (i(0), x(0), r(0),p(0)) = 0, with p(0) ^ 0. 
Note that for fixed t(0), x(0) and p(0), we have m distinct choices for r(0), equal to the m distinct 
real roots of er m (t(0), x(Q), T,p(0)). These choices for r(0) give m distinct waves that travel in 
different directions. The curve (t(s), x(s)) in physical space is the space-time ray that Gaussian 
beams are concentrated near. For a proof that the Gaussian beam construction is only possible 
near this ray, we refer the reader to |30j . 

The following lemma summarizes some results related to the Hamiltonian flow above. We will 
use the second point to argue that changing variables s t is always allowed. The last point is 
needed in the proof of the non-squeezing lemma (Lemma |3|). 

Lemma 1 Let (t(s),x(s),T(s),p(s)) be a null bicharacteristic of the Hamiltonian flow (0) with 
initial data such that \p(0)\ 7^ and a~ m (i(0), x(0), r(0),p(0)) = 0. Without loss of generality, 
assume that the parametrization is taken so that t(0) > 0. Then for s € [0, 00), we have 

1- o- m (s) = a m (t(s),x(s),T(s),p(s)) = 0, 

2. t(s) — t(0) is strictly increasing and 

t( s )-t(o)>c 1 r ' . v m = 1 ' 

W W " \log (1 + CaWOjp- 1 *) , m>l, 
where C\ and C2 are independent of s and initial data. 

3. There is a constant A independent of s and initial data such that 

\p{s)\ > |p(0)|e- A ( t ^- t (°)). (10) 
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Proof: We compute 



. da m - da 

cr m (s) = —^-t + 



dt dy 
da m da m da. 



. , da n 
x H — - — r ■ 



da n 



8r 
da m 



dp 
dam da 



■ P 



dc 



da Ti 



= 0, 



dt dr dy dp dr dt dp dy 

and since a m (0) = 0, we have that a m (s) = 0, proving the first point. 

Since a(x,t,p,r) — 0, we note that |p(so)| = implies that t(sq) = 0. Thus, if |p(so)| = then 
p(so) — ^§^-(so) = 0. Hence, we have that for all s 6 [0, oo), \p{s)\ ^ by uniqueness for solutions 
of ODEs and |p(0)| 7^ 0. Recall that for a polynomial with distinct root, both the polynomial 
and its derivative cannot vanish at the same point. Since \p{s)\ 7^ and strict hyperbolicity 
imply that a m (t(s), x(s), r,p(s)), as a polynomial of r, has distinct roots and a rn = on the null 
bicharacteristic, we have that ^^(s) 7^ 0. Continuity implies that i(s) — ^r-(s) never changes 
sign and, by the choice in paramctrization, we have that t(s) > 0. Hence, t(s) ~ t(0) is strictly 
increasing for all s € [0, 00). 

We next show that there is a constant C such that |r(s)| < C|p(s)| for all s and initial data. 
This is obviously true if r = 0. When t 7^ 0, let £ = p/r and observe that, by homogeneity, 

= a m (t,y,T,p) = T m a m (t,y, 1,£) , 

on the null bicharacteristic. Hence £ is a root of the polynomial equation 

m 

1+ £ 9p(t,y)^ = 0. 

I/3| = 1 

Since the coefficients gp{t, y) are bounded it follows that there is a constant such that |£| > C > 0, 
which proves that \t(s)\ < C\p(s)\ for all s. 
We can now bound p(s) as follows. 



\p\ 



da 7 , 



dy 

m— 1 



m— 1 



E E H 

J=° \\P\=m-j J 
m— 1 

< C E E \p\ m \rY = cY^ E \p\ m \tr <a\p\ m , 

3=0 \/3\=m-j 3=0 \P\=m-j 

for some constant c\ independent of s and initial data. Let A ~ c\/cq where cq is the constant in 
©. Moreover, set t(s) = t(s) - t(0) > and note that t(s) = i(s) > 0. Then, 



£ W .)|V*» 



2 p ■ p e 



2\t 



2Xt\p\ 2 e 2M > -2\p\\p\e 



2x1 



lm+1 2At 



2Ac |pr +1 e 2At = 



> -2 Cl \p\'--e- 

Consequently, |p(s)| 2 e 2At(s) > |p(0)| 2 and (10]) follows. With A m = (m - 1)A, we obtain 



t = 



da m {t,y,T,p) 
dr 



> co\p\ m - 1 > cob(0)| 



It follows that i(s) > cqs for m = 1, since Ai = 0. For m > 1, 



scob(0)| 



m— 1 



< 



t( S ')e AmtV W 



t~0) 



d0 = 



,A m t(s) 



m— 1^ — A TT1 i(s) 

£( 8 ') e *»*(«'W 
- 1 
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This proves the stated logarithmic growth of t(s) for m > 1. □ 

As mentioned above, an immediate consequence of this lemma is that we can use t to parametrize 
the Hamiltonian flow © instead of s, since the lemma guarantees that for a fixed to £ [0, T), there 
exists a unique sq such that t(so) = to- With a slight abuse of notation we will now write x(t) and 
p(t) for the Hamiltonian flow parametrized by t. 

Following [30] . we define a phase function <f> and amplitude functions a 3 via Taylor polynomials. 
After changing variables s — > t in the formulation used in [30] we can write: 

fc+i , , fc+i , 

<Kt,v) = E ^Mt)y p = Mt) + yp(t) + y-M(t) y + J2 ^M*)/ - 

|/3|=0 P ' " l/3|=3 P ' 

fc-2j-l 

%•(*,*/) = E «T a ^(*)/- (ii) 

|/3|=0 

We can now define the preliminary /c-th order Gaussian beam Vk(t,y) as: 

5 fc (t, tf )= E e i a J -(*,i/-i(*)) e < *( t '»-»'W^. 

Applying the operator P to this beam and collecting terms containing the same power of e, we 
have 

./ \ 



Pv k (t,y)= V e r c r (t,y))e^y-^^ , 



where c r (t,y) are smooth functions independent of e. The construction in [30] then proceeds to 
make c r (t, y) vanishe up to order k — 2(r+rn) + f on y = x(t). To obtain this, (x(t),p(t)) must follow 
the Hamiltonian flow and the coefficients in the Taylor polynomials should satisfy ODEs, which 
are given as follows. By assumption |(H2)| above, whenever p ^ 0, we can define m Hamiltonians 
Hi(t,x,p) implicitly by the relations cr m (t, x, —Ht(t,x,p),p) — for £ — 0, ...,m — 1. For any 
choice H — Hi, the first several ODEs are 

x = d p H(t,x,p) , p = -d y H(t,x,p) , (12) 

<t> = -H + p ■ d p H(t, x, p) , M = —A — MB - B T M - MCM , 

with 

_ d 2 H _ d 2 H c _ d 2 H 

dy 2 ' dpdy ' dp 2 

By (|TU| the Hamitonian will be well-defined for all times if p(0) ^ 0. Moreover, we have the 
following result for the Hessian matrix M(t) of the phase that guarantees that the leading order 
shape of the beam stays Gaussian for all time: 

Lemma 2 (Ralston '82, |30j) Suppose that M(0) is chosen so that it has a positive definite 
imaginary part, Af(0)±(0) = p(0), and M(0) = M T (0). Then, for all t £ [0,T], M(t) will be such 
that M(t) will have a positive definite imaginary part and M(t) = M T (t). 

Before we can fully define a Gaussian beam, one last point needs to be addressed. Since x(t) 
and p{t) are real, if 0o(O) is real and M(0) chosen as in Lemma[5] then the imaginary part of 4>(t, y) 
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will be a positive quadratic plus higher order terms about x(t). Thus, we must only construct the 
Gaussian beam in a domain on which the quadratic part is dominant. To this end, we use a cutoff 
function p v g C°°(IR n ;]R) with cutoff radius < r\ < oo satisfying, 



1 for \z\ < r\ 



for < r] < oo 



,<„(: ).: and ,,„[ :) = { ° f ° r N " 2?? . (13) 

1 for i] = oo 

Now, by choosing ry > sufficiently small, we can ensure that on the support of p n (y — x(t)), 
^s(j>{t, y — x(t)) > S\y — x{t)\ 2 for t € [0, T]. However, note that for first order beams the imaginary 
part of the phase is quadratic with no higher order terms, so that the cutoff is unnecessary. Thus, 
to include this case we let the cutoff function be defined for 77 = 00 by p^ = 1. We are now ready 
to finally define the fc-th order Gaussian beam Vk(t,y) as: 

rfe/21-1 

v k (t,y)= 4(»-«(*))«i(*.I/-«(*))« i * M)/e . (14) 

3=0 



2.2 Superpositions of Gaussian Beams 

In the previous section, we introduced Gaussian beam solutions that satisfy Pu — in an asymp- 
totic sense (see [30] for a precise statement), without too much concern for the values of the solution 
at t = 0. The initial value problem that a single Gaussian beam Vk{t,y) approximates has initial 
data that is simply given by its values (and the values of its time derivatives) at t = 0. While they 
resemble the initial conditions for ([2]), they are quite different since for example the phase, </>, for 
Vk is complex valued and Vk is concentrated in y. 

Our goal is to create an asymptotically valid solution to the PDE ©, thus we must also 
consider the m distinct pieces of initial data given in the form of time derivatives of u at t = 0. To 
generate solutions based on Gaussian beams that approximate the initial data for @, we exploit 
the linearity properties of P. That is, we use the fact that a linear combination of two Gaussian 
beams, with different initial parameters, will also be an asymptotic solution to Pu = 0, since each 
Gaussian beam is itself an asymptotic solution. Building on this idea, we take a family of Gaussian 
beams that is indexed by a parameter z £ Kq, where Kq is the compact subset of R" discussed 
in the introduction that contains the support of the amplitudes of the initial data for ([2]). We 
will use the notation xe (t; z), pi(t; z), <fie(t, y — x(t; z); z), etc, to denote the dependence of these 
quantities on the indexing parameter z and on the m choices for H{t,x,p) = Hf (t,x,p) denoted 
by £ = 0, . . . , m — 1. Thus, we will write the fc-th order Gaussian beam Vk,i(t, y; z). Note that the 
cutoff radius 7/ may vary with z and £ between beams, however, as the beam superposition is taken 
over the compact set Kq and < I < m — 1 , there is a minimum value for r\ that will work for all 
beams in the superposition. Thus, we form the fc-th order superposition solution Uk(t,y) as 

m_1 / 1 \? f 

where the phase and amplitude that define the Gaussian beam, 
r*/21-i 

v klt (t, y-z)= ]T eip v (y - x e (t; z))a td {t, y - x e (t; z); z ) e ^t,v-Mt^)/e , 
3=0 
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are given by 

1 k+1 1 

<j>t(t, V, z ) = 0£ : o(i; z)+y pi(t; z) + y -Mg(t; z)y + ^ —^{t] z)y fs , (16) 

\/3\=3 P - 

k-2j-l 

aej(t,y;z) = ^ —a iJtP (t;z)y^. (17) 

|/3|=0 P - 

We remind the reader that each Vk^it, y \ z) requires initial values for the ray and all of the amplitude 
and phase Taylor coefficients. The appropriate choice of these initial values will make Uk(0,y) 
asymptotically converge the initial conditions in $2$. The first step is to choose the origin of the 
rays and the initial coefficients of (f>£ up to order k + 1. Letting the ray begin at a point z and 
expanding $(y) in a Taylor series about this point, 

fe+i . 

*(f)= J^e( z )(y ~z) + error , 

101=0 1 

we initialize the ray and Gaussian beam phase associated with Tg as 
x e (0;z) = z, Pe(Q;z) = V^O) , 

A&(0;«) = 3j$(z)+ild„xn, fo/KO; = |j8| = , |/3| = 3, . . . , k + 1 . 

Determining the initial coefficients for the amplitudes involves a more complicated procedure. As 
with the phase, we expand all of the amplitude functions in Taylor series, 

fe-2j-l 

^j(y) = XI o|^j>/3( z )(y ~ 2 )' 3 + error ■ 

101=0 

Next, we look at the time derivatives of Uk at t = 0, to the lowest order in e. We equate the 
coefficients of (y — z) to the the corresponding terms in the Taylor expansions of Bg j and recalling 
that 

d t (f>e(0,0;z) = <j>£ fi (t;z) -xz(0;z) -pe{0;z) = -H<t(0,Xi(0;z),pe(0;z)) =: t £ , 
we obtain the following m x m system of linear equations: 



1 


1 




a OA0 








(jT m _l) £ 










(iTo)™- 1 • 


• (iT m -l) m - l _ 




_Om-l,O,0_ 




_Am-l,0,p_ 



Since the Tg are distinct, this Vandermonde matrix is invertible, so the solution will give the initial 
coefficients for the first amplitude for each of the m Gaussian beams. If we proceed with the next 
orders in e, we would obtain the same m x m linear system for [dcy^g, . . . , a,m—i,j,p] T > except that 
the right hand side will not only depend on the Taylor coefficients Bgj^, but also on previously 
computed af, 9 , 7 , q < j, coefficients and their time derivatives. Thus, all of the necessary initial 
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coefficients for each of the m Gaussian beams can be computed sequentially. Summarizing this 
construction, we have that at t = and I = 0, . . . , m — 1, 

8tu h (0,y) = (Jj^j 2 e- e J^sip v (y-z)b^(y)e l ^/z dz + o(e°°) , (18) 

where <j>{y) is the Taylor expansion of &(y) + i\y — z\ 2 /2 to order k + 1 and each bg t j is the same 
as the Taylor expansion of Bgj up to order k — 2j — 1. The 0(e°°) term is present because some 
of the time derivatives fall on the cutoff function p^y — xi(t; z)). The contributions of such terms 
decays exponentially as e — > 0, since the derivatives of prj(y — xe(t; z)) are compactly supported 
and vanish near y = z. 

Remark 2.1 For ease of notation and exposition, in 0) we have taken the phase $ to be the same 
for all of the m initial data, however, this is not a requirement. We can form m different Gaussian 
beam superpositions that each satisfy one of these m conditions with a specific phase and the rest 
with zero initial data. Then, summing these m solutions we obtain a more general solution of @) 
with m different phase functions for each of initial data piece. 

Remark 2.2 In the initialization of M^(0; z) we take its imaginary part to be given by ild nxn 
for simplicity. All of the results in this paper can be carried out if we instead took the imaginary 
part to be given by i"/ld nxn , for some constant 7 > 7 and adjusted the normalization constant in 
\15)) appropriately. However, it is important to note that the constants throughout this paper will 
depend on 7 and that in general, we expect that increasing 7 will increase the evolution error in 
the Gaussian beams and that decreasing 7 will increase the error in approximating the initial data. 

This completes the construction of the Gaussian beam superpositions Uk for the initial value 
problem (|2|) for general m-th order strictly hyperbolic operators. From this point, we will assume 
that the parameter n is chosen as n = 00 for the first order superposition u\ and that for higher order 
superpositions it is taken small enough (and independent of z and £) to make 3^(t, y— xi (t; z); z) > 
S\y — Xi{t;z)\ 2 for t 6 [0,T] and y € supp{p v (y — Xg(t;z))}. Lemma S] below shows that this is 
always possible. 



2.3 The Schrodinger Equation 

The construction of Gaussian beams for hyperbolic equations can be extended to the Schrodinger 
equation by replacing the operator P with a semiclassical operator P £ . Then, we can similarly 
construct asymptotic solutions to P 6 u — as e — > 0. In this section, we briefly review the 
construction presented in [33] for the Schrodinger equation (0} with a smooth external potential 
V(y). Note that the small parameter e represents the fast space and time scale introduced in the 
equation, as well as the typical wavelength of oscillations of the initial data. 
We recall that the fc-th order Gaussian beam solutions are of the form 

rk/21-1 

v k (t,y;z) = J2 ^p^y-x^z^^y-x^z^z)^^-^'^, 
j=o 

where the phase and amplitudes are given in (ITC1) and (JT7J) and p v is the cutoff function (p~5|) . 
Furthermore, the subindex has been suppressed since for the Schrodinger equation there is 
only one choice for H{t,x,p), namely 



H(t,x, P ) =*L-+V{x) 



(19) 
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The system (TT21 for the bicharacteristics (x(t; z),p(t; z)) is then given by 
x = p , x(0; z) = z , 

P=-VyV, P(0]Z) = Vy$(z) . 

The equations for the phase and amplitude Taylor coefficients are derived in the same way as for 
the strictly hyperbolic equations. The phase coefficients along the bicharacteristic curve satisfy 

\p\ 2 

fo = ^--V{x{t)) , 
M = -M 2 - d 2 y V{x{t)) , 

<(>? = - £ (7 _i )!(/3 _ 7) , ^^-7 + 2 " dfr |/3| = 3, . . . , k + 1 . 

The amplitude coefficients are obtained by recursively solving transport equations for aj^ with 
< k — 2j — 1, starting from 

ao.o = --^a 0fi Tr(M(t; z)). 
These equations when equipped with the following initial data 

<M0;z) = $o(z) , M{0;z) = d 2 <S>(z) + iU nxn , 

fo{0;z) = $p(z), \p\ = 3,...,fe + l , aj , p {0;z)=A jt p(z) , 

have global in time solution, thus Vk(t, y; z) is well-defined for all < t < T. 
The fc-th order Gaussian beam superposition is finally formed as 



Uk(t,y) = ( - — ) / v k (t,y;z)dz . 



\2ire J 



Ko 



As in the case of strictly hyperbolic PDEs, we will assume that the cutoff parameter r\ is chosen 
as T) = oo for the first order superposition u\ and that for higher order superpositions it is taken 
small enough (and independent of z) to make 9</>(£, y — x(t; z); z) > S\y — x(t; z)\ 2 for t £ [0, T] and 
y € supp{p^(y — x(t; z))}. Again, Lemma 0] ensures that this can be done. Furthermore, we note 
that Remark 12.21 concerning the initial choice of the imaginary part of M(0; z) also applies to the 
superposition for the Schrodinger equation. 



3 Error Estimates for Gaussian Beams 

In this section we prove the asymptotic convergence results for superpositions of Gaussian beams 
given in our main result, Theorem[T] The corner stone of our error estimates are the well-posedness 
estimates for each PDE. Since they are crucial to our analysis we summarize them here. 

Theorem 2 The generic well-posedness estimate 

\\u(t, < KO, Oils + Ce« [ \\G[u](t, 0|| ia dr , (20) 

Jo 



applies to 
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• the m-th order strictly hyperbolic PDE JJ|) with = P, q = and \\-\\s the Sobolev space-time 
norm, 

m—l 

l| 9 tVM|| H m-*-l , 

£=0 

where H s is the Sobolev s-norm (H = L 2 ), 

• the wave equation f3[) with Q — d 2 — c(y) 2 A, q = 1 and || • \\s the e-scaled energy norm, 

IKv)l|£ "(Ti^F + |Vw|2dy ) 1/2 ' (21) 

• and the Schrddinger equation with = P £ , q = — 1 and || • ||s £/ie standard L 2 norm. 

Proof: The results for the wave and Schrodinger equation are standard and can be found in most 
books on PDEs. The result for m-th order equations is a bit more technical to prove and appears 
in Section 23.2 of [I3J (Lemma 23.2.1). □ 

Remark 3.1 Since the wave equation is a second order strictly hyperbolic PDE, we have two 
distinct well-posedness estimates in terms of two different norms. Furthermore, we note that \\-\\e 
is only a norm over the class of functions that tend to zero at infinity, which we are considering 
here. 

When Theorem [5] is applied to the difference between the Gaussian beam superposition, u^, 
and the true solution, u, for any one of the PDEs that we are considering, we obtain the following 
estimate for t £ [0,T], 

H(t, •) - u(t, 0||s < IK(0, •) - u(0, Oils + Ce* f ||9[u fc ](T, Oil , (22) 

Jo 

with the appropriate choices for 0, q and || • \\s- We will refer to the first term on the right hand 
side as the error in approximating the initial data or the initial data error and to the second term 
as the evolution error. 

Using the ideas in [34] , we prove Theorem [4] which shows the convergence rate in e of the 
Gaussian beam superposition to the initial data for any given Sobolev norm. Thus, this theorem 
extends a result of |34j , so that we can use it to estimate the initial data error in the more general 
well-posedness estimates above. 

The evolution error has been estimated in the work of previous authors [35] [53] 03 [32] . The 
necessary steps used by those authors are quite general and can be applied to any strictly hyperbolic 
equation as well as any linear dispersive wave equation as long as it is semi-classically rescaled. 
Following these ideas, we show in Lemma [5] that for all of the PDEs that we consider, Q[uk] can 
be written in the form 

,/ 

©K] = ^'^^(fla,*,,/^^) +0(£°°) , Tj > 0, fi e L 2 {K Q ), (23) 
3=1 

so that the key to estimating the evolution error is precise norm estimates in terms of e of the 
oscillatory integral operators Q aj , gj , v ■ L 2 (K ) M> L°°([0, T]; L 2 (R ra )), defined as follows. For a 
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fixed t £ [0,T], a multi-index a, a compact set Kq C E", a cutoff function ,o r) (fTB")) with cutoff 
radius < 77 < 00 and a function g(t,y;z), we let 

(Q a , g , v w)(t,y) (24) 

:= e- 2 ^ f w(z)g(t, y; z){y - x(t; z)) a e^-^^ p n (y - x(t; z))dz , 

with the functions g(t, y; z), <f>(t, y; z) and x(t; z) satisfying for all t £ [0, T], 

(Al) x{t;z)€C°°([0 ) T\xK ), 

(A2) <j>{t,y;z),g(t,y-z) £ C°°([0,T] x E" x 

(A3) V (/>(£, 0; z) is real and there is a constant C such that for all z, z' £ Ko, 

|V v 0(t, 0; z) - V v 0(t, 0; z')| + |ar(i; z) - x(t; z')\ > C\z - z'\ , 

(A4) for \y\ < 2rj (or for all y if 77 = 00), there exists a constant 5 such that for all z £ Kq, 

^(t,y;z)>S\y\ 2 , 

(A5) for any multi-index /3, there exists a constant such that 

sup |c^g(t,y;z)| < C/3 . 
zeK 

With this definition, the following norm estimate of Q a , g ,ri will be proved in Sectional 
Theorem 3 Under the assumptions (Al)-(A5), 

SUp \\Q a .g. V \\ L 2 < C(T). 

se[o,T] 

This theorem improves the norm estimate of Q a ,g.r) given in |221 123] . which has an additional 
factor e f in the right hand side, with 7 = (n — l)/4 for the wave equation and 7 = u/4 for 
the Schrodinger equation. Instead of estimating the integral directly, we follow the arguments in 
[31] to relate the estimate of the oscillatory integral to the operator norm, through the use of 
an adjoint operator. An essential ingredient in estimating the operator norm is the non-squeezing 
lemma (Lemma 13]), which states that the distance between two physical points is comparable to 
the distance between their Hamiltonian trajectories measured in phase space, even in the presence 
of caustics. Using Theorem [3] we are able to prove the same convergence rate for Gaussian beam 
superpositions over physical space that is achieved in [5] for beam superposition carried out in full 
phase space. Thus, we improve on the error estimates given in (22j [23] to obtain error estimates 
for the Gaussian beam superposition that are independent of dimension as given in Theorem [Q 
We conclude this section with two remarks. 

Remark 3.2 The assumption of C°° smoothness for all functions is made for simplicity to avoid 
a too technical discussion about precise regularity requirements. In this sense, Theorem [7] and 
Theorem^ can be sharpened, since they will be true also for less regular functions. 

Remark 3.3 // the condition in assumption (A4) is satisfied for all y there is no need for the 
cutoff function in the definition of the operator Q in \2J$ . We treat this case by taking r\ = 00 and 
defining = 1. The operators with r\ — 00 are used in the case of first order Gaussian beams. 
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3.1 Gaussian Beam Phase 

In this section, we show that the Gaussian beam phase <f> given in (1161) is an admissible phase for 
the operators Q a , g ,r)- We begin with a lemma based on the regularity of the Hamiltonian flow map, 
St, stating that the difference \z — z'\ is comparable to the sum \p(t; z) — p(t; z') \ + \x(t; z) — x(t; z')\. 
Note, however, that because of caustics it is not true that \z — z'\ is related in this way to either 
of the individual terms \p(t; z) — p(t; z')\ or \x(t; z) — x(t; z')\. 

Lemma 3 (Non-squeezing lemma) Let St be a Hamiltonian flow map, 

(x(t; z),p(t; z)) = S t (x(0; z),p(0; z)), 

associated to a strictly hyperbolic PDO US\) or to the Schrodinger operator \19i) . Let Kq be a 
compact subset of R n and assume that p(0; z) is Lipschitz continuous in z € Kq for the flow 
associated with the Schrodinger operator. Additionally, assume that inf zS if |p(0;z)| = 6 > for 
the flow associated with the strictly hyperbolic PDO. Under these conditions, there exist positive 
constants C\ and C2 depending on T and 5, such that 

ci\z - z'\ < \p(t; z) - p(t; z')\ + \x(t; z) - x(t; z')\ < c 2 \z - z'\ , (25) 

for all z, z' G Kq and t G [0, T]. 

Proof: We prove the result for the flows associated with the two types of operators separately, 
however, we use the common notation, 

Z = (z, Po ) = (z,p(0; z)) , Z 1 = (z',p' Q ) = (z',p(0; z 1 )) , 

X = S t [Z) = (x,p) = (x(t;z),p(t;z)) , X' = S t {Z') = {x',p') = {x(t; z'),p(t; z')) . 

We begin with the flow associated with the Hamiltonian for the Schrodinger operator. Let us 
introduce the set K.q = {(z,p(0;z)) : z G K } and note that S t is invertible with inverse S-t and 
regular for all t so that 

sup sup 

te[o,T] zeconv(K; ) 

where conv(i?) denotes the convex hull of the set E. Now, noting that 

X - X' = £ ±S t {sZ + (1 - s)Z') d s = f Q ^{ S Z + {l-s)Z') [z _ z>)ds (26) 

and taking l\ norms, since sZ + (1 — s)Z' G conv(i^o)> we have 

\x - x'\ + \p -p'\ = \\X - X'\\x < C\\Z - Z% = C(\z -z'\ + \p - p' \) < C'\z - z'\ , 

where we have used the Lipschitz continuity of po in z. This gives the right half of (|25p . By the 
equivalent of for S'_ t , we have 

\z ~ A <\z- z'\ + bo -Pol = \\Z- Z'\\i <C\\X- X'\\i = C(\x - x'\ + \p-p'\) , 

which completes the proof of the lemma for the Hamiltonian flow associated with the Schrodinger 
operator. 



dSt(Z) 
dZ 



< C 



sup sup 

t£[0,T] X£conv(S t (K )) 



dS- t (X) 
dX 
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For the flow associated with the strictly hyperbolic operator, we follow the same idea, but we 
have to be careful near \p\ = 0. Thus, in addition to /Co, we introduce the sets 

B s = {{z,p) : z€K , \p\ < 6}, 

JCq = conv (/Co U Bs) ■ 

Note that St is regular away from |po| =0 by Lemma Q] as \p(t)\ 0, for all t so 

sup sup 

te[o,T] zeK \B s/2 

Thus, again by (|2T>1) we obtain, 

\x-x'\ + \p-p'\ = \\X -X'^ <C\\Z - Z'\\ 1 = C{\z - z'\ + \p - p' \) <C'\z - z'\ , (27) 

provided that p(s) = (1 — s)po + sp' satisfies info< s <i \p{s)\ > 5/2, which guarantees that sZ + (1 — 
s)Z' £ IC \ B S / 2 for < s < 1. On the other hand, suppose that inf < s <i \p(s)\ < 6/2. We define 
p* as the point on the line connecting p to p' with smallest norm and let s* = argmin 0<s<1 \p(s)\ 
so that p* = p(s*). Then 

->\p*\ = \ Po - s*(p Q -p' Q )\ > \p Q \ - s*\p Q — Pol > - bo -Pol- 
Hence, |p -p' | > S/2. Now, let 

d = sup sup ||St(Z)||i < oo, 
te[o,T]ZeK 

so that 

HX-X'Hi < ||X|| 1 + ||X / || 1 <2d<24bo-Pol<C|l^-^lli ! 

o 

which shows that (J27J) holds for all Z, Z' £ JC - This gives the right half of (Jl5j). For the left half 
of ([25]) , we note that by Lemma [lj in£te[Q,T],zeK \p(t': z )\ > <5 > 0, so that we can consider the 
inverse map S-t in exactly the same way as St above and show that \\Z — Z'\\\ < C\\X — -X"'||i- 
Thus we obtain the left half of (j25|) . 

\z - z'\ <\z- z'\ + |po - Pol = \\Z- Z'Wt < C\\X - X% < C(\x - x'\ + \p-p'\) . 

Thus, the proof of the lemma is complete. □ 

We are now ready to show that the phase function for a fc-th order Gaussian beam is an 
admissible phase for the operators Q ajffj?? . 

Lemma 4 Let <p(t, y; z) and x(t; z) be the phase and central ray associated to a k-th order Gaussian 
beam for fH), f3J) or ^ as given in Section ] 2.2\ and Section \2.S\ The rays x(t; z) satisfy (Al) and 
<p(t,y;z) satisfies assumptions (A2) through (A4) iff] is sufficiently small. In the case of k = 1, rj 
can take any value in (0, oo]. 

Proof: The smoothness assumptions (Al) and (A2) follow from smoothness of initial data and 
smoothness of the coefficients in the underlying PDE. By definition, V y 4>(t, 0; z) = p(t; z) and (A3) 
follows from the non-squeezing lemma, Lemma |3] Finally, since the lower order terms of <fi are real, 

fe+i . 

^(t,y;z)=y(QM(t;z))y+ £ -^(t; z)y? . 

101=3 P ' 



dSt(Z) 



dZ 



< C. 
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Recalling that by Lemma [U 3M(i; z) is positive definite, we therefore have that when \y\ < 2 77, 

k+l 



mt,y;z)>C\y\ 2 - £ i||3<^(t; 



W\ 



101=3 

>C\y\ 2 -\y\ 2 J2 ill^(*;«)IU-(2T7)^l- a 

|/3|=3 P " 

> |y| 2 [C-277 £ i||^/s(<^)|U-(2»7)l fl l- 3 ) 



>^)|y| 2 , 

where the constant S(rj) is positive for small enough 77 and independent of t G [0, T], since ^ (t; z) 
are smooth functions of t. This shows that <p(t,y;z) satisfies (A4). When k = 1, there are only 
quadratic terms in the phase and in this case 4>(t, y; z) will satisfy (A4) for any choice of 77 € (0, 00]. 
□ 

3.2 Representations of P[uk] in Terms of Q^g^ 

In this section, we show that several of the intermediate quantities in the proof of Theorem [T] can 
be written as sums involving the operators Q Q)ffir? . 

Lemma 5 Let P be an m-th order strictly hyperbolic operator and P £ a semiclassical Schrddinger 
operator, satsifying the assumptions stated for (dp ( in Section \2. 1\) and ^ respectively. Let Uk be 
the corresponding Gaussian beam superpositions given in Section \2.2\ and Section \2.3[ Then P[uk] 
and P e [uk] can be expressed as a finite sum of the operators Q: 

1 x § J ( lj > k/2 - m + 1 

Z7r S j=i { l 3 > fc/2 + 1 

with xk {z) the characteristic function on K$ and Q ajl g, iri satisfying assumptions (A1)-(A5) if 77 
is sufficiently small. In the case of k = 1, 77 can take any value in (0, 00]. 

Proof: For a single Gaussian beam Vk,i{t,y; z) for the hyperbolic operator, following the 
discussion in |3D] and Section l!Q1 we have 

\k/2\-l 

P[v k , e (t, y; z)} = £r P^V - x ^ *))<V(*> V\ z)e^*™^)^/ £ + E k (t, y; z) , 

r— — m 

where E k contains terms that are multiplied by derivatives of the cutoff function. For first order 
beams 77 = 00 and poo = 1 making E\(t,y;z) = 0. For higher order beams, E k = in an 77 
neighborhood of xg(t;z) and since 77 is small enough so that the imaginary part of <f>n is strictly 
positive for 77 < \y — xi(t; z)\ < 277, E k will decay exponentially as e — > 0. Thus, Ek — C(e°°) for 
all orders of beams and t G [0, T]. 

Each c r can be expressed in terms of the symbols of P: 

c- m +ji{t,y, z) = Laej-i + a m (t,y,d t 4>e,^ y 4>e) a>e,j + Rt,j(t,y, z) , 
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where aij = for j ^ [0, \k/2\ — 1] and L is given by (using the Einstein summation), 



T . ( da m - da 

V 9 Pj d Vj 



with V = (<9t,V), y — (t,y) and p = (r,p). The function Re t o(t, y; z) = and the functions 
Re,j(t, y; z) for j > are complicated functions of the m — 2 and lower order symbols of P and the 
functions <f)£, ag.fi, • • ■ , and their derivatives. We note that since agj(t,y; z) are compactly 

supported in z € Aoj so are the functions c r j(t, y; z). 

By the construction of a Gaussian beam, c,-^ vanishes up to order k — 2(r + m) + 1 on y = 
(t; z). Note that fc — 2(r + m) + 1 may be negative, in which case, c r ,g is not necessarily on 
y = xg(t; z). Thus, by Taylor's remainder formula, 

Cr,t(t>V> z ) = E Cr,£, a (t,y;z) (y - x e (t;z)) a , 

|a| = [fc-2(r+m)+2] + 

for some coefficient functions c r ,e, a {t,y; z) that are compactly supported in z € Kq and [a]+ = 
max[a, 0]. For the superposition Uk(t,y), we have 



P[u k ] 



2tT£ J 



^ / P[v k ^(t,y;z)}dz 
ftm-l /f*/2l-l 



£=0 \ r——m 



EE/ e r /^(!/-^(*;«))^(* > l/;«)e i *' (t *-* <(t; * }; * )/e ds , 



modulo additive terms that are 0(£°°). Now substituting for c r .^ we can rewrite this expression in 
terms of the operators Q, 

/ 1 m-1 / ffc/21-1 ' 



V27re 



E E 



E 



2tt 



£=0 \ r=-m \ |a| = [fc-2(r+m)+2]+ ' 

x c. Aa (t, y; z)(y - a*(t; z ))« e <*i(t,»-i(t;»);*)/e dz 

-1 / ffc/21-l / 



E e r+H/2 (Q a , CrJ ,^XK )(t,y)\\ , 

£=0 \ r=-m V |a|=[fc-2(r+m)+2] + // 

modulo additive terms that are 0(£°°) and where XAoM ^ s ^ ne characteristic function on K and 
the functions Cr t g, a (t, y; z) satisfy condition (A5) since the coefficients of P, (f>e,p and Q>i,jf} are all 
smooth. By Lemma|4]the operators also satisfy (A1)-(A4) under the condition on r\. 

To simplify the notation, let j = 1, . . . , J < oo enumerate all of the combination of i, r, and a 
in the triple sum above and rewrite the sums as 

J 



P[u k ] = 



with lj > (fc/2 — m + 1) and gj = Cr-g-a-. Thus, we have the desired result for P. 

Now, for each Gaussian beam Vk(t,y;z) for the Schrodinger equation defined in Section |2.3 
following [53], we compute 

"ffc/21-i 



+ Ek(t,y;z) 
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where again Ef. = 0(e°°). Note that 



e -* E F[ae* e ] =aG{t,y)-ieLa- —Aa , 



with 



V(y) and L = d t + V y <j) ■ V y + -A y(l 



Thus, we have 



P' 



\k/2]-l 

E 



e j ai e i,p/£ 



rfe/21-i 

^ ' e- 7 a,jG — ieLdj — — A y aj 
j'=o 



J=0 

where for convenience aj = for j £ [0, |7c/2] — 1] and 

d r = a r G — iLa r -i — -A y a r -2 , T = 0, . . . , \k/2] + 1 . 

By construction of the phase coefficients in Section [2731 we have that on y — x(t; z), G vanishes to 
order k + 1 and by construction of the amplitude coefficients, the quantity —iLa r ^\ — ^ A y a r _2 
vanishes to order k — 2r + 1. Thus, d r vanishes to order k — 2r + 1 on y — x(t; z), where again 
we remark that if k — 2r + 1 is negative, d r does not vanish on y = x(t; z). By Taylor's remainder 
formula, we have 

d r (t,y;z) = ^2 dr, a (t,y;z) (y - x{t;z)) a , 

|a| = [fc+2-2r] + 

with d rt0l (t, y; z) compactly supported in z G Kq. Hence, 



|7c/2]+l 

p £ [vk] = Pv E £T 

r=0 

For the superposition, we have 



£ d r<a {y-x{t;z)Y\e i ^ + 0{e°°) 

. |a| = [fc+2-2r] + 



P e [u k ] 



V2vre 



'ifo 
ffe/21-1 

£ £ £ r+|tt|/2 (Q Q ^, Q ,,XKo)(^2/) + 0(£ 00 ) • 

r=0 |a| = [fc+2-2r] + 



As above, the smoothness of V(y), <fip and a,- p guarantees that d r ^ a satisfies (A5), while LemmaU 
ensures that also (A1)-(A4) are true for Q a ,d r , a ,r) un der the condition on r]. Again, we let j = 
1, . . . , J < oo enumerate all of the combination of r and a in the double sum above. Thus, 



= ( ^ ) E ^' (2<™ **o) (*, 2/) + 0(O 



with > (fc/2 + 1) and gj — d rj aj . Thus, the proof is complete. □ 
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3.3 Initial Data Errors 

In addition to establishing the role of the Q-operators and estimating their norms, we need to 
consider the convergence of the Gaussian beam superposition to the initial data to be able to 
prove Theorem [TJ We start with the following lemma. 

Lemma 6 Let <f> € C°°(Ko) be a real-valued function and Aj £ Cq°(Ko). With k > 1, define 

N 

N 



uk(y) 



where cf>(jj — z; z) is the k + 1 order Taylor series of $ about z, dj(y — z\ z) is the same as the 
Taylor series of Aj about z up to order k — 2j — 1 (but may differ in higher order terms) and p^ 
is the cutoff function \13\) with < i] < oo . Then for some constant C s , 



\Uk - u\\ H , 



Proof: The proof of this lemma is based on a discussion in [33]. We first assume that r\ < oo. 
Looking at each of the terms in the sum above separately, the estimate depends on how well 

e 3 A J {y)e^ {y)/e =(—Y [ e j AMe^'^-^ ' 2e dz (28) 

V 27r£ / ^R" 



is approximated by 

^ 2 J e 3 Pn (y - z) aj (y - z; z y<t>iv-^)/s-\ y -zf '/2e dz ^ 



Since we want to estimate the difference between these two functions in a Sobolev norm, we need 
to consider differences between their d@ derivatives in the L 2 norm. Since derivatives of (|29[) 
that fall on the cutoff function p n (y — z) vanish in a neighborhood of z = y and the integrand is 
compactly supported in y and z, they will be 0(e°°) in the L 2 norm and do not contribute to the 
estimate. Thus, when we differentiate ||2SJ) and O by d$ and pair the results from the sequence 
of differentiations, the terms that will contribute to the estimate will be of the form 

(db) 2 ^c^dyiMy)} ( II d v m [my)-\y-*m) (so) 



for (gS} and 



( — 



ei-tcirtgh [aj{y _ z - z)] JT d ^ -z;z)-\y- z\ 2 /2] (31) 



for ([2^)l . where C^ 7 ^ are combinatorial coefficients. These terms are integrated over R" in z and 
summed over the multi-indexes Si +7 = /3, the index I = 1 . . . I7I, and multi-indexes |7i|, . . . , \^t\ > 
1, 71 + . . . + 7^ = 7. Furthermore, we have that max m [7 m ] < [7] — I + 1. These formulas can 
be obtained through long but straightforward calculations. The important part is to recall that 
icf>(y — 2; z) — \y— z\ 2 /2 is the k + 1 order Taylor series of i$(y) — \y — z\ 2 /2 about z and aj(y — z; z) 
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agrees with the Taylor series of Aj(y) about z up to order k — 2j — 1, so that (13"TT) will agree with 
the Taylor series of ([50)1 up to order 



k — 2j — 1 — k + 1 — max |7„ 

Km<£ 



>fc-2j-l-|5i|-| 7 |+i 

= fc-2j-l-| J 8|+l. 
Thus, we have that \\d!?(uk — u)\\l 2 can be estimated by a sum of terms of the form 



{toe) 



i<S>(y)/e~\y-z\ 2 /2e 



Be(y)e 



Pv (y - z)b e (y - z; z y*(y-^)/e-\y-A 2 1^ 



dz 



(32) 



where I < \0\, \B e - b e \ = 0{\y - z |k-2:H0l+<) and |$ - (f>\ = 0{\y - z\ k+2 ). Now, the proof of 
Theorem 2.1 in [34] can be applied directly to (1321) to obtain the estimate, 



N 



3=0 



Thus, we have the result for r\ < oo. The extension for p^ = 1 follows directly, since the cutoff p v 
for v < oo introduces 0(£°°) errors in the L 2 norm: 



[1 - p n (y - z)] b e (y - z; z)^^-^'^^ ' 2e dz 



0(e c 



L- 



as 1 — p v vanishes in a neighborhood z = y and the integrand is compactly supported in z. □ 

Using Lemma [6l we can estimate the asymptotic convergence rate of the superposition solution 
to the initial data. 



Theorem 4 For the Gaussian beam superposition, Uk, given in Section \2.2\ and the solution, u, 
to the strictly hyperbolic PDE (0), we have 



\d e t u k (0,-)-d e t u(0,-)\\ Hs <Ct 



,e 2 



-e-s 



for some constant C# )S and < £ < m — 1. 

Similarly, for the Gaussian beam superposition, u k , given in Section [2~B and the solution, u, to 
the wave equation {3J) ; we have 

\\u k (0,-)-u(0,-)\\ E <Cei . 



Furthermore, for the superposition, Uf., given in Section \2.S\ and the solution, u, to the Schrddinger 
equation Q), we have 

||« fc (0,-)-«(0,-)|| La <Ce* . 

Proof: The proof of this theorem for hyperbolic PDEs follows directly from Lemma since 
for each power of e, dfiik{0,y) and dfu(0,y) given in (fl8|) and ([2]), respectively, are exactly in the 
assumed form in the Lemma [6] 

The result for the wave equation follows after noting that 

|K(0, •) - u(0, Oils < Ce (|K(0, •) - u(0, -)\\ H i + \\dtu k (0, ■) - d t u(0, -)M . 

Similarly, the result for the Schrodinger equation follows directly from the definition of the Uk and 
u at t = in Section 12.31 □ 
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3.4 Proof of Theorem Q] 

We prove the results for each type of PDE separately. For the strictly hyperbolic m-th order PDE 
@, applying the well-posedness estimate given in Theorem [5] to the difference between the true 
solution u and the fc-th order Gaussian beam superposition, u k , defined in Section |2.2[ we obtain 
forte[0,T], 

m—1 

Pt ') _ U *(*> Olllffm-*-! 

1=0 

(771 — 1 'J 1 \ 

J2 ||af[«(O ) .)-u fc (O,0]|| flW _ < _ 1 + / ||PM(r,.)|| L2 dr . 
fco - 70 / 

The first term of the right hand side, which represents the difference in the initial data, can 
be estimated by Theorem |4] and the second term, which represents the evolution error, can be 
estimated by Lemma [5] to obtain 

m— 1 

J2 1 1 °t [«(*.■)- «* (*>•)] 1 1 



«=0 



<C(T) £ #-™+i+^^ sup ||Q Q , )9i ^|| i2 



j=l *6[0. T 1 



with £j > (fc/2 — to + 1) and Q ajj g jjr) satisfying (A1)-(A5), for small enough rj when A; > 1. Thus, 
using Theorem [3l we obtain 



m — 1 



£ ||#[u(v)-u*(v)]||^-*-i 



which completes the proof for strictly hyperbolic PDEs. Since the wave equation is a second order 
strictly hyperbolic PDE, applying the above estimate to ([3]), we obtain for t G [0, T], 

l 

\Ht r )-u k (t r )\\ E <E^\\d^[u{t,-)-u k (t,-)]\\ H1 _ t <C{T)ei , 

1=0 

which completes the proof of Theorem [1] for the wave equation. 

For the Schrodinger equation ((4]) , applying the well-posedness estimate given in Theorem [2] to 
the difference between the true solution u and the fc-th order Gaussian beam superposition, Ufc, 
defined in Section [2~3l we obtain for t € [0,T], 

\\u k {t, •) - U(t, -)\\ L 2 < \\u k (0, •) - u(0, 0IU= + " / \\P"[u k }(T, -)\\ L 2dT . 

£ Jo 

The initial data part of the right hand side can be estimated by Theorem [4] to obtain ||u(0, •) — 
u k(0, ")IU 2 < Ce^ . With the help of Lemma we can estimate the second part of the right hand 
side as 



T J 
- / Hi^ttfcKT.OIMT^yV'- 1 SUp || Q c 

£ Jo =1 te[o,T] 



,9j 



J L2 + 0(e°°) , 
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with £j > (fc/2 + 1) and, as above, Q ajt g j!V satisfying (A1)-(A5), for small enough r\ when k > 1. 
Again, using Theorem [3] and combining, we obtain, 

,/ 

\\u k (t,-)-u(t,-)\\ L 2 <Cs% +e*X;<7 a (T) + 0( e ~)<:<7(T>4 . 

i=i 

Thus, the proof of Theorem Q] is complete. 

4 Norm Estimates of Q a ,g,r/ 

In this section we prove Theorem [3] We follow the ideas in [SJ to relate the estimate of the 
oscillatory integral to the operator norm, through the use of an adjoint operator. A key ingredient 
in estimating the operator norm is the non-squeezing lemma (Lemma [3]), which allows us to obtain 
a dimensionally independent estimates for the oscillatory integral operator. 

4.1 Operator Norm Estimates of Q a , g ,ri 

We let Q* g „ be the adjoint operator and consider the squared expression, 

{Q*a,g,r,Q<x,g, V u)(t, Z) 

= e -n-\a\ f u ( z !y l 4>(t,y~x(t-,zyy)/e el4 , { t,y- x (t-,z):z)/s g (^ y . y . ^ 

JU n xK 

x (y - x(t; z')) a {y - x(t; z)) a p v (y - x(t; z'))p n (y - x(t; z))dydz' 
:=e- n - lal f I £ a Jt,z,z')u(z')dz', (33) 



9^ 

K 



where 

Il g {t, z, z') = / e W*^(t*)J)/e e W,v-<M*)rt/e g {t, y - z')g(t,y;z) 

x (y - x(t] z')) a (y - x(t; z)) a p v (y - x(t; z')) Pri {y - x(t; z))dy 
el Mt, v ,z,z')/e g ^ y+ -. z ' )g ^ y + i ;z) 

x (y - &x) a (y + Ax) a p v (y - Ax) p v (y + Ax) dy , 
after a change of variables and 

x(t; z) + x(t; z') 



x — x(t, z, z) := 
Ax = Ax(t, z, z) := 



2 

x(t; z) — x(t; z') 



">P(t,y,z,z') := (j>(t, y + Ax; z) - <j>{t\ y - Ax; z) 

This symmetrization will simplify expressions later on. 
Recall Schur's lemma: 

Lemma 7 (Schur) For integrable kernels K(x,y), 

J K(x,y)u(x)dx < ^sup J \K(x,y)\dy^j ^sup J \K(x,y)\d 
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Using Schur's lemma, we can now deduce that 



I Qa,g,v I II 2 — SU P 

w£L 2 (K ) 



w\ 



L- 



< sup 

weL 2 (K ) 



\Q* a ,g,r,Q<*,9,r,M\L* 
IMU a 



< e 



< e 



-I«l 

-\a\ 



sup / |ij Jt,z,z')\dz') sup / |/£ (*,;>:, z')|eb 



z'SKo Jk 



sup / |^. g (t,z,z')|dz' 



upon noting that |I£ )ff (t, z, z')\ = \I e ag {t,z\z)\. 
Before continuing, we need some utility results. 



4.1.1 Utility results 

We will prove a few general results that will be useful in the proof of Theorem [3J 

Lemma 8 (Phase estimate) Let rj be the same as in assumption (A4). Then, under the as- 
sumptions (A2)-(A4), t G [0, T], and y such that \y ± Ax\ < 1r\ (or all y if rj — oo ), we have: 

• For all z,z' S Kq, there exists a constant 5 independent oft such that 

^(t,y,z,z')> U [\y + Ax\ 2 + \y - Ax\ 2 ] = S\y\ 2 + ±8\x(t; z) - x(t; z')\ 2 . 

• For \x(t; z) - x(t; z')\ < 9\z - z'\, 

inf \V y TP(t,y,z,z')\>C(6,n)\z-z'\ , 

where fl(t,fi) — {y : \y — Ax\ < 2/i and \y + Ax\ < 2(i} and C(9,fi) is independent oft and 
positive if 9 and fi < r\ are sufficiently small. 

Proof: By assumption (A4), there exists a constant 5 independent of t such that 

3V (*, V, z, z 1 ) = 3f0(t, y + Ax; z 1 ) + 3$(f, y - Ax; z) > S (\y + Ax\ 2 + \y - Ax\ 2 ) 

= 2S\y\ 2 + \s\x-x'\ 2 . 



For convenience, we divide by 1/2 to eliminate the factor in front of 6\y\ 
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For the second result, we have 

\V y iP(t,y,z,z')\ > \$tV y iP(t,y,z,z')\ 

= |HV B 0(t, y + Ax; z') - HV„0(t, y - Ax; z)\ 

= m y 4>{t,Q;z')-m y <j){t,Q;z) 

+ m y <f>{t, y; z') - 0; z') - y; z) - HV w 0(t, 0; «)] 

+ y + Ax; z') - y; z 1 ) 

- [RV„0(t, y - Ax; z) - y; z)] 

> |5ftV y 0(t, 0; z') - m y (f>{t, 0; z)| 

- | (V„0(t, y; z') - V w 0(t, 0; z')) - (V„0(t, y; z) - V v ^(t, 0; z))\ 

- |V w 0(t, y + Ax; z') - V w 0(t, y; z')| 

- |V y 0(t, y - Ax; z) - V„0(t, y; z)\ 
—\ E\ — E2 — E^ — E^ . 

Using assumption (A3) we have for E\, 

E x = |SV s #,0;z')-SV^(t,0;2)| = |V B ^(t, 0; z') - V„0(t, 0; z)| 
> C|z- z'| - |x- x'|, 

where x = x(t; z) and x' = x(t; z'). To estimate E 2 , we first note that on Sl(t, fi), 



1 1 A 1 1 A 
—y Ax H — y H — Ax 

2 y 2 2 tf 2 



<-(|y-Ax| + |y + Ax|)<2 A1 . 



Then, by the Fundamental Theorem of Calculus and the smoothness assumption (A2), for y e 
Q(t, ji) we have, 



£-2 



/ [d^sy^^-d^sy^^yds 
Jo 



<C\z-z'\\y\<CMz-z'\ , 



with C\ independent of t € [0,T]. Similarly, with y <G 0(t, /i) and C2 independent of t € [0,T] and 
s G #0, 

|V„4>(i,y± Ax;s) - V y 0(t,y;s)| < C 2 |Ax| = hj 2 \x-x'\ , 

which shows that E% + E% < C 2 \x — x'\. Using these estimates for the case \x — x'\ < 0\z — z'\ we 
then obtain 



|V W V(*. V, z > z ')\ > C\z - z'| -\x- x'| - dii\z - z'| - C 2 |x - x'| 
> C\z - z'| - (1 + C 2 )0\z - z'| - Cm\z - z'| 
=: C(6,n)\z-z'\ , 

where C(0,/z) is independent of t £ [0,T] and positive if and /U are small enough. □ 
Next we have a version of the non-stationary phase lemma. 
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Lemma 9 (Non-stationary phase lemma) Suppose that u(y; Q € C^°(D x Z), where D and 
Z are compact sets and ip{y\C,) € C°°(0) for some open neighborhood of D x Z . If V y ip never 
vanishes in O, then for any K = 0, 1, . . 



\a\<K 

where Ck is a constant independent of 



<C K e K V ( 1^(^01 e^^ )/£ dv 



Proof: This is a classical result. A proof can be obtained by modifying the proof of Lemma 7.7.1 
of [H]. However, we omit the details for the sake of brevity. □ 
With this lemma in hand we can estimate \I^ g (t, z, z')\. 

Lemma 10 Under the assumptions (A2)-(A5), for any K = 0, 1, . . fixed < p < r\ < oo, s > 
and t €E [0,T], there are constants Ck and C s independent oft such that 

\I e a Jt,z,z')\ <C K e n ' 2+ ^—— ± ' 7m +C s e s , (34) 

9 1 + mf^o^^) \Vyip{t, y, z, 0/Ve| K 

where Q(t, fi) = {y : \y — Ax\ < 2p and \y + Ax\ < 2p} C {y : \y\ < 2/i} is a compact set. 

Proof: By the definition of I e a (t, z, z'), we have 

F a Jt, z, z') = f e^^'^git, y + x; z')g(t,y + x;z) 

x (y - Ax) a (y + Ax) a p n (y + Ax)p v (y - Ax)dy 

e i^t,y,z,z'ye g ^ y+ -. z i )g{tjy + s . ;z ) 



x(y- Ax) a (y + Ax) a p n {y + Ax)p n {y - Ax)dy 

e i^t, v ,z,z')/e g ^ y + -. Z ' )g{ti y + S; . z) 

n(t,ri)\n(t,fj,) 

x(y- Ax) a (y + Ax) a p v (y + Ax)p v {y - Ax)dy 

The integral I\ will correspond to the first part of the right hand side of the estimate in the lemma 
and I2 to the second part. We begin estimating 1%. By Lemma [8] and (A5), for a fixed t, we 
compute, 

< C / \y- Ax\^\y + Ax\^e- s ^ Ax ^y +Ax ^/ E dy . 

Now, using the estimate s p e~ as2 < (p/e) p / 2 a~ p / 2 e~ a!>2 ^ 2 , with p = \a\, a = S/e and s = \y — Ax\ 
or \y + Ax\, and continuing the estimate of 1%, we have for a constant, C, independent of t, z and 

\h\<c( £ -) M ( e -A(i y+ A,i 2 +i,-A,n dy 

< C ( £ -) lal [ e -l\y\ 2 -i\^\ 2 d y<Ce n / 2+ ^e-^ 2 . 
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Thus, we have proved the needed estimate for I\ for the case K = as well as the the case K > 
when hrfygntt./i) |VyV(^j 2/? z ? z ')l = 0- Therefore, in the remainder of the proof we will consider the 
case K ^ and miy^u,^) | V J/ ?/'(i, y, z, z')| ^ 0. In this case, Lemma [5] can be applied to I\ with 
C = (t, z, z') G [0, T] x x to give, 



l/J|<# 
L8 <K \ 



l<9|/2 



x / |fl£ [{y - Ax) a {y + AxYg'gp+p-] \ e~^dy 

\0\<K V ' ' ' \/3i+/3 2 = ( 9"' s H*'W 
/3i<2a 

where — p v (y ± Ax), u(t,z,z') ~ infj, e n(t.^) | VyV'C^j J/i z > Z ')/V^\ an d Cr- is independent of t, 
z and z' . By assumption (A5) and since p^ is uniformly smooth and t, z, z' vary in a compact 
set, \d@ 2 [g'gp^p^] | can be bounded by a constant independent of y, t, z and z' . We estimate the 
other term as follows, 

\d^[(y-Ax) a (y + Ax) a }\ <C £ | (y - Ax) a ~^ (y + Ax) a ~^ \ 

/3n+/3i2=/3i 

<C \V~ Ax| |a| " l/3111 \y + Ax| |a| - |tel . 



/3n+/3i2=/3i 
/3n,/3i2<a 



Now, using the same argument as for the K = case, we have 

\d& [(y - Ax) a (y + Ax) a }\ \d$* [g'gp+p-] \ e~^dy 



0(t,^i) 



<C I \V- Ax| |Q| - |/Jl11 \y + Ax\W~^e- C ^ /E dy 

(9n,/3i2<a 

< C(fe)£ " +l " | - | ^J + l " | - | ^ l e -^|A a; r = C . (/32)e n/2 + |a|-|^ 1 |/2 e -| |A^ 



and consequently, 



l/3|/ 2 

— /2+|a|-|/3 1 |/2 p -f|Ax| ;i 

i/(t,3,*') Mf-|/ ' 1 

|/3|<-fsT V ' ' ' h+h=P 
Pt<2a 



< CK£ „/2+|a| e -||Ax|= £ 1 



Kt,2,2') 2iC " l/31 
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Using the fact that will be bounded by the minimum of the K — and K > estimates, we 
have 



17x1 < C£ n / 2 +H e -tl A *l 2 mm 



W\<K 



1 



z/(t,z,z') 2K ~ l/3 ' 



Noting that for positive a, 6, and c, 

min[a, 6 + c] < min[a, b] + min[a, c] and min[l, 1/a] < 2/(1 + a) , 



we have, 



l/3|<^ 



;/(i,z,z') 2K -^l 



< min 

l/9|<* 



1- 



1 



< 



St 

l/3|<-ff 



i/(t,z,z') 2Ar " l/31 
2 



< C K - 



+ v{t, z, z') 2K -\P\ ~ * 1 + i/(t, -2) *') K 



This shows the I\ contribution to the estimate (fM|) . It remains to show the smallness of I2 ■ Indeed, 
since either \y + Ax| > /i or \y — Aa;| >/i on 17(i, 77) \ f2(t, /i), we get in the same way as for I\ in 
the K = case, 



U2I < c 



(I) 



e 2i 



-(ly+A^ + ly-AxI 2 ) 



Sl(t,jj)\Sl(t,p) 



for any s > with C s independent of i. This concludes the proof of the lemma. □ 



4.1.2 Proof of Theorem U 

We now have all of the ingredients to complete the proof of Theorem [3] We fix t E [0, T] and start 
with the estimate, 



\Q a ^\\h < e- n - H (sup / \r a Jt,z,z')\dz') , 



K 

derived in the beginning of Section 14.11 and we turn our attention to estimating the integral of 
\I e a<g (t,z,z')\. 

We will use the shorthand notation Ify(t,z,z') = Xd(z, z')I e a g (t, z, z'), where xd(z,z') is the 
characteristic function on a domain D C Kq x Kq. We will consider two disjoint subsets of Ko x Ko 
given by 

Dx{t, 6) = {(z, z') : \x(t; z) - x(t; z')\ > 0\z - z'\} 
D 2 (t,0)={(z,z') : \x(t;z)~x{t;z , )\ < 9\z - z'\} , 

where 9 is the small parameter 9 in Lemma HJ Note that D\ U D2 = Ko x Kq and D\ n D2 = 0. 
Thus, 



\r a jt,z,z')\dz' 



K 



\r Di (t,z,z')\dz' + / \r D2 {t,z,z')\dz' . 

K JKo 
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The set D\ corresponds to the non-caustic region of the solution. There, we estimate If) (t, z, z') 
by taking K = and s = n + \ a\ in Lemma ITUl 

\r Dl (t,z,z')\dz' < Ce n ' 2+ ^ [ e~^ x ^- x (t;z ' )|2 +e n/2 dz' 

K 1 J K 



< C£ n/2+\ a \ f e -€\*-*'\ 2 dz' + C\K \s n+ \ a 
JKo 

< Ce n/2+\ a \ / s n-l e -^ ds + Ce n+\a\ 

Jo 



< Ce n+H . 

The set corresponds to the region near caustics of the solution. On D2, we estimate If, (t, z, z') 
using Lemma [TO] again where we pick \i small enough to allow us to use also Lemma [5] Letting 
R = sup 2 Z 't=K \ z — z'\ < 00 be the diameter of Kq and s — n + \ct\, we compute 



n . r ( -i;\x(t:z)-x(t:z')\ 2 \ 

Ih 2 (t,z,z') dz' <C £ -+H / — — — — ——- +e *)dz! 



<Cr e i+H f 1 + Ce n +I Q I 

«o 1 + 



C(9,m)|z-z'I X K 



<C^+I Q I / — . . — s^ds + Cg"+l a l 

< Ce n+|a| , 

if we take K = n + 1. 

Since all of the constants are independent of the fixed i € [0, T], by putting all of these estimate 
together we obtain 

\\Q a , g J\h <£~"- |Q| ( sup / <C, 
for all t £ [0, T], which proves the theorem. 



5 Numerical study of convergence 

In this section, we perform numerical convergence analyses to study the sharpness of the theoretical 
estimates in this paper for the constant coefficient wave equation with sound speed c(y) = 1, for 
which H(t,x,p) = ±|p|. The ODEs that define the Gaussian beams are solved numerically using 
an explicit Runge-Kutta (4,5) method (MATLAB's ode45). We use the fast Fourier transform 
to obtain the "exact solution" and use it to determine the error in the Gaussian beam solution. 
When the norms require it, we compute derivatives via analytical forms rather than numerical 
differentiation. 



5.1 Single Gaussian Beams 

First, we study the convergence rate for a single Gaussian beam to show the sharpness of the 
estimate proved in [30]. For 2D, this estimate states that for a single Gaussian beam Vk(t,y), 



\\(d?-c(y) 2 A)v k (t,y)\\ L 2<Ce k ^-^ , 
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for t £ [0, T]. Using the well-posedness estimate for the wave equation, we obtain the rescaled 
energy norm estimate 

\\v k (t, •) - u(t, -)\\ E < |M0, •) - u(0, -)\\ E + T e sup \\{d 2 t - c(y) 2 A)v k (t, , 

S£[0,T] 

for t G [0, T], where u is the exact solution to the wave equation. Taking the initial conditions for 
u to be the same as the 1-st, 2-nd and 3-rd order Gaussian beams at t = (modulo the cutoff 
function), we obtain the asymptotic error estimate, 

IMV)-u(V)IU <Ce fc /2+i/2 _ 

To test the sharpness of this estimate, we investigate the numerical convergence as follows. We let 
the Gaussian beam parameters for the 1-st order Gaussian beam be given by 



x(0) = 
M(0) = 






i 
2 + i 



p(0) = 
ao,o(0) = 1 



Mo) = o 



and pick H = — \p\. The additional coefficients necessary for the Taylor polynomials of phase and 
amplitudes for 2-nd and 3-rd order Gaussian beams are all initially taken to be 0. Note that even 
though the phase and amplitudes for 1-st, 2-nd and 3-rd order beams are the same at t = 0, their 
time derivatives will not be, as the ODEs for the higher order coefficients are inhomogeneous. 
Thus, since the initial data for the exact solution, u, matches the initial data for the Gaussian 
beam, u depends on the order of the beam. With this choice of parameters, we generate the 1-st, 
2-nd and 3-rd order Gaussian beam solution at t = {0.5, 1}. For 2-nd and 3-rd order beams we 
additionally use a cutoff function with rj = 1/10. At each t, we compute the rescaled energy norm 
of the difference v k — u. The asymptotic convergence as e —> is shown in Figure [T] We draw the 
attention of the reader to the following features of the plots in Figure [U 

1. For 1-st order beam: The error decays as e . 

2. For 2-nd and 3-rd order beams: The error for larger e is dominated by the error induced by 
the cutoff function. In this region, the error decays exponentially fast. As e gets smaller, the 
error decays as e 3 / 2 for the 2-nd order beam and e 2 for the 3-rd order beam. 

The numerical results agree with the estimate given in [30 and, thus, the estimate for single 
Gaussian beams is sharp. 

5.2 Cusp Caustic 

We consider an example in 2D that develops a cusp caustic. The initial data for u at t — is given 

by 

u(0,y) = e -i°MV(-^i)/- . 
Thus, the initial phase and amplitudes are given by 

Hv) = -Vx+Vl , A , (y) = e- 10 ^ 2 , A OA (y) = 0. 
For the initial data for u t (0,y), we take 

u t (o,y) = f l -My) [AoAy) + eAoAy)] + [Ao,oAv) + £4>,m0/)]) , 
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Figure 1: Single Gaussian beams: Asymptotic behavior of Vk — u in the rescaled energy norm 
at t = {0.5,1} for k = {1,2,3} order beams. The results are shown on log-log plots along with 
cie 1 , C2£ 3 / 2 , and c^e 2 to help with the interpretation of the asymptotic behavior. The asymptotic 
behavior agrees with the analytical estimates: 1-st order beam is 0(e 1 ), 2-nd order beam is 0(e 3 ^ 2 ), 
and 3-rd order beam is 0{e 2 ). 



where the $ t , A 0j o,t and A Q ^ it are obtained from the Gaussian beam ODEs and various y derivatives 
of <&, A 0i o and A ^. Specifically, we take $t(y) = +|V y $(j/)| so that waves propagate in the positive 
2/1 direction. As was shown in [34], this particular example develops a cusp caustic at t = 0.5 and 
two fold caustics for t > 0.5. 

To form the Gaussian beam superposition solutions, it is enough to consider Gaussian beams 
governed by H = — \p\, because of the initial data choice. We take the initial Taylor coefficients to 
be 

x(0;z)= M , 

<Po (0; z) = -z x + z\ , 

^(0;z) = 0, |0|=3,4, 
oi,o = , 

where only the necessary parameters are used for each of the 1-st, 2-nd and 3-rd order beams. For 
2-nd and 3-rd order beams, we additionally use a cutoff function with 77 = 1/10. We propagate the 
Gaussian beam solutions to t = {0,0.25,0.5,0.75,1}. At each t, we compute the rescaled energy 
norm of the difference Uk — u. The asymptotic convergence as e — > is shown in Figure [2] 
We draw the attention of the reader to the following features in the plots: 

1. At all t: 

(a) The asymptotic behavior of 1-st and 2-nd order solutions is the same and of order 0(e 1 ). 

(b) The asymptotic behavior of 3-rd order solution is of order 0(e 2 ). 

(c) For large value of e the error for 2-nd and 3-rd order solutions is dominated by the error 
induced by the cutoff function. In this region the error decays exponentially. 



p(0;z) 



-1 

2z 2 



i 
2 + i 



M(0; z) 

aoA0^) = ^Ao(z), \/3\ =0,1,2 
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Figure 2: Superpositions of Gaussian beams for cusp caustic: Asymptotic behavior of it/. — u 
in the rescaled energy norm at t = {0,0.25,0.5,0.75,1} for k = {1,2,3} order Gaussian beam 
superpositions. The results are shown on log-log plots along with c\e l and C2£ 2 to help with the 
interpretation of the asymptotic behavior. 



(d) For midrange values of e, 2-nd order solutions experience a fortuitous error cancellation. 
This is due to the cutoff function. Changing the cutoff radius r\ shifts this region. 

2. At t — 0.5 and t > 0.5, the asymptotic behavior of the error is unaffected by the cusp and 
fold caustics, respectively. 

For this example we note that the convergence rate for odd beams k £ {1,3} is in fact half 
an order better than our theoretical estimates. This improvement was also observed and proved 
for a simplified setting in [57] , where the analysis shows that the gain is due to error cancellations 
between adjacent beams; it is therefore not present for single Gaussian beams. The result in [27] 
is for the Helmholtz equation and only concerns the pointwise error away from caustics. Here the 
numerical results indicate that the same improvement appears for the wave equation in the energy 
norm even when caustics are present. We conjecture that this gain of convergence is due to error 
cancellations as well and that it is present for all Gaussian beam superpositions with beams of odd 
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order k. Moreover, we conjecture that the theoretical result is sharp for beams of even order k, 
giving us an optimal error estimate in the energy norm of ©(e^/ 2 ! ) for all k. 

6 Concluding Remarks 

Gaussian beams are asymptotically valid high frequency solutions to strictly hyperbolic PDEs 
concentrated on a single curve through the physical domain. They can also be constructed for the 
Schrodinger equation. Superpositions of Gaussian beams provide a powerful tool to generate more 
general high frequency solutions. In this work, we establish error estimates of the Gaussian beam 
superposition for all strictly hyperbolic PDEs and the Schrodinger equation. Our study gives the 
surprising conclusion that even if the superposition is done over physical space, the error is still 
independent of the number of dimension and of the presence of caustics. Thus, we improve upon 
earlier results by Liu and Ralston 22, 25] , 
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